In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['figure.figsize'] = [10,7]
How to project 3D data to 2D maps¶
Projecting 3D data to 2D cannot be done without distortions
You can conserve either but never all:
- equal-angle/conformal: projected circles will stay circular
- equal-area: projected circles have same area
- equidistant: projected circles have same radii in conserving direction
In [2]:
ax = plt.axes(projection=ccrs.Orthographic())
ax.set_global()
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[2]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e5453e20>
Cylindrical projections¶
- longitudes and latitudes are parallel
- longitudinal distances are equally-spaced and conserved
Angle conserving: Mercartor projection¶
In [3]:
ax = plt.axes(projection=ccrs.Mercator())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[3]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6e57d90>
Area conserving: LambertCylindrical¶
In [4]:
ax = plt.axes(projection=ccrs.LambertCylindrical())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[4]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6daefa0>
Distance conserving: PlateCarree projection¶
In [5]:
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[5]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d84340>
Good compromise without conservation: Miller projection¶
In [6]:
ax = plt.axes(projection=ccrs.Miller())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[6]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d512b0>
Pseudocylindrical/Elliptical projections¶
- typically are area conserving (but not all!)
- do not have parallel latitudes
- can have parallel longitudes
Area and distance conserving: Sinusoidal projection¶
In [7]:
ax = plt.axes(projection=ccrs.Sinusoidal())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[7]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d185e0>
Area conserving: Mollweide projection¶
- longitudes are parallel and unequally-spaced
- meridians are ellipses
In [8]:
ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[8]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6c869a0>
Compromise with less distortion: Aitoff projection¶
- not equal-area, based on azimuthal equidistant projection
- longitudes are not parallel for less distortion at the edges
- edge is 2:1 ellipse
In [9]:
ax = plt.axes(projection=ccrs.Aitoff())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[9]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6c51d00>
Area conserving with less distortion: Hammer projection¶
- equal-area, based on azimuthal equal-area projection
- longitudes are not parallel for less distortion at the edges
- edge is 2:1 ellipse
In [10]:
ax = plt.axes(projection=ccrs.Hammer())
ax.stock_img()
ax.gridlines()
ax.tissot()
Out[10]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e4427730>
Conclusion¶
Use area conserving elliptical projections like Mollweide or Hammer
In [11]:
# All you need to create projections is contained in the coordinate reference system
import cartopy.crs as ccrs
# In order to plot a map, one has to tell matplotlib which projection to use
ax = plt.axes(projection=ccrs.Mollweide())
# This is not an Axes object but a GeoAxes, a subclass of Axes representing a projection
print(type(ax))
# All normal matplotlib plotting methods work but have to include a transform in which format the data was defined
data = np.arange(4).reshape([2,2])
ax.imshow(data, transform=ccrs.PlateCarree())
<class 'cartopy.mpl.geoaxes.GeoAxes'>
Out[11]:
<matplotlib.image.AxesImage at 0x7fb8e438e580>
In [12]:
# Be aware that cartopy uses the Earth ellipsoid as bases for projections
# You can redefine the globe to be a perfect sphere
ax = plt.axes(projection=ccrs.Mollweide(globe=ccrs.Globe(ellipse='sphere')))
Arrays can be plotted normally with imshow¶
In [13]:
# Define projection
ax = plt.axes(projection=ccrs.Mollweide())
# Define grid where first index is latitudes and second index is longitudes
data = np.arange(72).reshape([12,6])
# As we defined a PlateCarree grid, we have to use the transform accordingly
# Imshow fits the grid assuming the data is regularly spaced from -180° to 180° longitude and -90° to 90° latitude
# However, imshow always flips the y-axis (origin='lower')
imap = ax.imshow(data, transform=ccrs.PlateCarree())
# Include colourbar
plt.colorbar(imap, location='bottom')
Out[13]:
<matplotlib.colorbar.Colorbar at 0x7fb8e4306700>
Creating regular interpolated data of a spherical shell¶
In [14]:
# Load in grid data
# In this case this is 3d grid data on a regular cartesian grid
# However interpolation should work regardless of shape and size
density = np.load('rho.npz')['arr_0']
print(f'Density array has shape {density.shape}')
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_xlabel('x [cells]')
ax1.set_ylabel('z [cells]')
ax1.imshow(np.log10(np.nansum(density, axis=1).T), origin='lower')
ax2.set_xlabel('x [cells]')
ax2.set_ylabel('y [cells]')
ax2.imshow(np.log10(np.nansum(density, axis=2)).T, origin='lower')
plt.tight_layout()
Density array has shape (159, 159, 159)
In [15]:
# Index grid
ind = np.indices(density.shape)
# Set origin to centre
centre = (np.array(density.shape)-1)/2.0
ind = (ind.T - centre).T
# Calculate radius
radius = np.sqrt(ind[0]**2 + ind[1]**2 + ind[2]**2)
# Cut out radial shell at given radius, here +-1 cell is enough -> flatten arrays
radcut = 75
mask = (radius <= radcut+1) & (radius >= radcut-1)
rad = radius[mask]
rho = density[mask]
x = ind[0][mask]
y = ind[1][mask]
z = ind[2][mask]
# Calculate longitude in degree
lon = np.rad2deg(np.arctan2(y, x))
# Calculate latitude (90 degree north) in degree
xyrad = np.sqrt(x**2 + y**2)
lat = np.rad2deg(np.arctan2(z, xyrad))
# Define equally-spaced interpolation points
step = 1.0
alon = np.arange(-180+step, 180, step)
alat = np.arange(-90+step, 90, step)
# Create grid from points
ilon, ilat, irad = np.meshgrid(alon, alat, radcut)
# Interpolate data
from scipy.interpolate import griddata
# Griddata interpolates the data with given points on the left to the new points on the right
irho = griddata((lon, lat, rad), rho, (ilon, ilat, irad))
In [16]:
# Define projection
ax = plt.axes(projection=ccrs.Mollweide())
# Plot interpolated data
# Centre is 0° longitude and 0° latitude
# Origin='lower' as imshow starts at top
# Cartopy is defined for geospatial data where the east is to the right
# In astrophysics we want to look from the 'inside-out' where east is to the west
# Therefore, we have to flip the longitude with extent
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))
# Include colourbar
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')
Out[16]:
<matplotlib.colorbar.Colorbar at 0x7fb8dd15d760>
In [17]:
# Interpolation does not know about cyclic data
ax = plt.axes(projection=ccrs.Mollweide())
# Plot shifted data, here centre is 180° longitude and 0° latitude
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(360,0,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')
Out[17]:
<matplotlib.colorbar.Colorbar at 0x7fb8dd06de20>
In [18]:
# Add cyclic data
large = np.where((lon>160) & (lon<180))
rho = np.append(rho, rho[large])
rad = np.append(rad, rad[large])
lat = np.append(lat, lat[large])
lon = np.append(lon, lon[large]-360)
small = np.where((lon<-160) & (lon>-180))
rho = np.append(rho, rho[small])
rad = np.append(rad, rad[small])
lat = np.append(lat, lat[small])
lon = np.append(lon, lon[small]+360)
# Redefine interpolation points as cyclic
step = 1.0
alon = np.arange(-180, 180, step)
alat = np.arange(-90+step, 90, step)
ilon, ilat, irad = np.meshgrid(alon, alat, radcut)
irho = griddata((lon, lat, rad), rho, (ilon, ilat, irad))
# Plot shifted data
ax = plt.axes(projection=ccrs.Mollweide())
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(360,0,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')
Out[18]:
<matplotlib.colorbar.Colorbar at 0x7fb8dcf5a250>
Plotting vector data¶
In [19]:
# Load velocity data
velazi = np.load('velazi.npz')['arr_0']
velpol = np.load('velpol.npz')['arr_0']
# Apply radial cut
vazi = velazi[mask]
vpol = velpol[mask]
# Add cyclic data
vazi = np.append(vazi, vazi[large])
vpol = np.append(vpol, vpol[large])
vazi = np.append(vazi, vazi[small])
vpol = np.append(vpol, vpol[small])
# Interpolate data
ivazi = griddata((lon, lat, rad), vazi, (ilon, ilat, irad))
ivpol = griddata((lon, lat, rad), vpol, (ilon, ilat, irad))
In [20]:
# Plot density map
ax = plt.axes(projection=ccrs.Mollweide())
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')
# Overplot quiver data
ax.quiver(-alon, alat, np.squeeze(ivazi), np.squeeze(ivpol), transform=ccrs.PlateCarree())
Out[20]:
<matplotlib.quiver.Quiver at 0x7fb8dcee9a60>
Use regrid_shape to regrid quiver data¶
In [21]:
# Regridding is a feature of cartopy which interpolates vector data onto regular points in projection space
# The input can be either an integer which defines the minimum number of vectors in one direction or directly given by a tuple of ints
# Plot density map
ax = plt.axes(projection=ccrs.Mollweide())
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')
# Overplot quiver data
ax.quiver(-alon, alat, np.squeeze(ivazi), np.squeeze(ivpol), transform=ccrs.PlateCarree(), regrid_shape=44)
Out[21]:
<matplotlib.quiver.Quiver at 0x7fb8dae7ce80>
Healpix¶
In [22]:
import healpy as hp
# Resolution is set by the parameter nside, which is generally a power of 2
NSIDE = 2
# Get number of pixels
print(f'The number of pixels is: {hp.nside2npix(NSIDE)}')
# Healpy assumes the pixels are stored in a certain order in a flattened array
data = np.arange(hp.nside2npix(NSIDE))
# Ring order assumes horizontal rings starting at the North Pole
hp.mollview(data, title="Ring order")
The number of pixels is: 48
In [23]:
# Nested order assumes horizontal rings starting at the North Pole
hp.mollview(data, nest=True, title="Nested order")
In [24]:
# Healpix comes with very powerful functions
# Standard coordinates are co-latitude and longitude
# Co-latitide is 0 at the North Pole, pi/2 at the equator and pi at the South Pole
# Longitude starts from centre, increases to the west and goes from 0 to 2pi
print(f'Co-latitude and longitude of pixel 0 in radians: {hp.pix2ang(NSIDE, 0)}')
# Get longitude and latitude of pixels in degree, Healpy goes from 0° to 360°
print(f'Longitude and latitude of pixel 0 in degree: {hp.pix2ang(NSIDE, 0, lonlat=True)}')
# Get unit vector of pixels
print(f'Unit vector of pixel 0: {hp.pix2vec(NSIDE, 0)}')
# Get resolution of nside in radians or arcmin
print(f'Distance between pixels at NSIDE={NSIDE} is about {np.round(hp.nside2resol(NSIDE, arcmin=True)/60, 2)} degree.')
# Get size of pixels in square radians or square degree
print(f'Size of pixels at NSIDE={NSIDE} is about {np.round(hp.nside2pixarea(NSIDE, degrees=True), 2)} square degree.')
# Get pixel numbers of overlap with different shapes
print(f'Indices of pixels overlapping with disc with radius 0.1 in radians at z-vector: {hp.query_disc(NSIDE, [0,0,1], 0.1, inclusive=True)}')
Co-latitude and longitude of pixel 0 in radians: (0.4111378623223478, 0.7853981633974483) Longitude and latitude of pixel 0 in degree: (45.0, 66.44353569089877) Unit vector of pixel 0: (0.2825970826302196, 0.28259708263021954, 0.9166666666666666) Distance between pixels at NSIDE=2 is about 29.32 degree. Size of pixels at NSIDE=2 is about 859.44 square degree. Indices of pixels overlapping with disc with radius 0.1 in radians at z-vector: [0 1 2 3]
Creating interpolated data of a spherical shell¶
In [25]:
# Increase the resolution to a similar one we used with cartopy
NSIDE = 64
# Define interpolation points, Healpy goes from 0° to 360°
ilon, ilat = hp.pix2ang(NSIDE, np.arange(hp.nside2npix(NSIDE)), lonlat=True)
# Shift data longitudes to be defined from 0° to 360°
slon = lon + 180
# Interpolate data
irho = griddata((slon, lat, rad), rho, (ilon, ilat, radcut))
# Plot data where we have to compensate for the 180° shift
hp.mollview(np.log10(irho), rot=(180,0,0), unit=r'density [g cm$^{-3}$]', title=None)
359.296875
Creating spherical integrations - how to distribute cell content¶
In [26]:
# Directly using cell centre gives bad results
# Index grid
ind = np.indices(density.shape)
# Set origin to centre
centre = (np.array(density.shape)-1)/2.0
ind = (ind.T - centre).T
# Calculate radius
radius = np.sqrt(ind[0]**2 + ind[1]**2 + ind[2]**2)
# Mask region
mask = (radius >= 1.0) & (radius<70.0)
rho = density[mask]
# Get pixel index of vectors
pix = hp.vec2pix(NSIDE, ind[0][mask], ind[1][mask], ind[2][mask])
# Sum up density in pixel
rhosum = np.zeros(hp.nside2npix(NSIDE))
rho = np.nan_to_num(rho)
np.add.at(rhosum, pix, rho)
# Plot result
hp.mollview(np.log10(rhosum), title='No distribution')
In [27]:
# Rough correction by spreading data to nearest neighbors, this does not work for all distances
# Sum up density in pixel
rhosum = np.zeros(hp.nside2npix(NSIDE))
rho = np.nan_to_num(rho)
np.add.at(rhosum, pix, rho/9)
# Add smoothed data by distributing cell content also onto next eight neighbors
pixneigh = hp.pixelfunc.get_all_neighbours(NSIDE,pix)
for i in range(8):
np.add.at(rhosum, pixneigh[i], rho/9)
# Plot result
hp.mollview(np.log10(rhosum), title='Simple smoothing onto eight neighbors')
In [28]:
# Better correction by getting the overlapping pixels with the cell
# Calculate angular extent of cells on the map in radians
# Here I estimate the mean angular extent of the regular cartesian grid to be sqrt(2) as we have all viewing angles
ext = np.arctan2(np.sqrt(2)/2, radius)
# Convert cell indices to vectors, flatten arrays
ext = ext[mask]
vec = np.array([ind[0][mask], ind[1][mask], ind[2][mask]]).T
# Get overlapping pixels and add data
rhosum = np.zeros(hp.nside2npix(NSIDE))
for i, ivec in enumerate(vec):
ioverlap = hp.query_disc(NSIDE, ivec, ext[i], inclusive=True)
np.add.at(rhosum, ioverlap, rho[i]/len(ioverlap))
# Plot result
hp.mollview(np.log10(rhosum), title='Simple pixel overlap')
In [29]:
# Further correction by weighting overlapping pixels with distance
# Calculate angular extent of cells on the map in radians
ext = np.arctan2(np.sqrt(2)/2, radius)
# Convert cell indices to vectors, flatten arrays
ext = ext[mask]
unitvec = (np.array([ind[0][mask], ind[1][mask], ind[2][mask]])/radius[mask]).T
ang = np.array(hp.vec2ang(vec, lonlat=True))
# Get overlapping pixels and add data
rhosum = np.zeros(hp.nside2npix(NSIDE))
for i, ivec in enumerate(unitvec):
# Get indices of overlapping pixels
ioverlap = hp.query_disc(NSIDE, ivec, ext[i], inclusive=True)
# Calculate distance of overlapping pixel centre to vector in radians
dist = np.arccos(np.dot(ivec, hp.pix2vec(NSIDE, ioverlap)))
# Calculate weight as some function of the inverse of distance
weight = (1/dist**0.4/np.sum(1/dist**0.4))
np.add.at(rhosum, ioverlap, rho[i]*weight)
# Plot result
hp.mollview(np.log10(rhosum), title='Weighted pixel overlap')